Electron acceleration due to high frequency 
instabilities at supernova remnant shocks 



M. E. Dieckmann,^ 1 K. G. McClements, b ' 2 S. C. Chapman, a 
R. O. Dendy, b ' a and L. O'C. Drury c 

O ; 

a Space and Astrophysics Group, Department of Physics, University of Warwick, 
^ ; Coventry, CV4 7AL, UK 

■ b EURATOM/UKAEA Fusion Association, Culham Science Centre, Abingdon, 

£ ■ Oxfordshire 0X14 3DB, UK 

j- J ■ Dublin Institute for Advanced Studies, 5 Merrion Square, Dublin 2, Ireland 

7-h ■ To appear in Astronomy and Astrophysics (accepted: February 15 2000) 

> ■ 

t>« | Abstract 

| Observations of synchrotron radiation across a wide range of wave- 

{Sj . lengths provide clear evidence that electrons are accelerated to rela- 

•O tivistic energies in supernova remnants (SNRs). However, a viable 

mechanism for the pre-acceleration of such electrons to mildly rela- 
tivistic energies has not yet been established. In this paper an electro- 
O ■ magnetic particle-in-cell (PIC) code is used to simulate acceleration 

^ \ of electrons from background energies to tens of keV at perpendic- 

$h ' ular collisionless shocks associated with SNRs. Free energy for elec- 

tron energization is provided by ions reflected from the shock front, 
| with speeds greater than the upstream electron thermal speed. The 

' PIC simulation results contain several new features, including: the 

^ ■ acceleration, rather than heating, of electrons via the Buneman insta- 

bility; the acceleration of electrons to speeds exceeding those of the 
shock-reflected ions producing the instability; and strong acceleration 
of electrons perpendicular to the magnetic field. Electron energization 
takes place through a variety of resonant and non-resonant processes, 
of which the strongest involves stochastic wave-particle interactions. 
In SNRs the diffusive shock process could then supply the final step 
required for the production of fully relativistic electrons. The mecha- 
nisms identified in this paper thus provide a possible solution to the 
electron pre-acceleration problem. 
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1 Introduction 



The detection of radio synchrotron emission from shell-type supernova rem- 
nants (SNRs) is a clear indication that electrons of typically GeV energy are 
being accelerated in such objects. There is now convincing evidence that syn- 
chrotron emission from some remnants extends to X-ray wavelengths (Pohl 
& Esposito 1998): this implies the presence of electrons with energies of order 
10 14 eV. The prime example of this is the remnant of SN1006, recent observa- 
tions of which using the ASCA (Koyama et al. 1995) and ROSAT (Willingale 
et al. 1996) spacecraft show that X-ray emission from the bright rim has 
a hard, approximately power-law spectrum. In contrast, emission from the 
centre is softer, with a strong atomic line component. The sharp edges and 
strong limb brightening observed at both X-ray and radio wavelengths in- 
dicate that: the acceleration site is the strong outer shock bounding the 
remnant; the acceleration is continuous; and the local diffusion coefficient of 
electrons near the shock front is substantially reduced relative to that in the 
general interstellar medium (Achterberg et al. 1994). The possibility that 
the X-ray emission from SN1006 is thermal bremsstrahlung has been exam- 
ined by Laming (1998), and found to be less tenable than the synchrotron 
interpretation. 

There is thus extensive observational evidence that the strong collisionless 
shocks bounding shell-type SNRs accelerate electrons to relativistic energies. 
The standard interpretation of extragalactic radio jet observations is also 
based on the premise that the relativistic electrons responsible for observed 
synchrotron emission are produced by shocks, although in this case the shock 
parameters are much less certain. Heliospheric shocks, on the other hand, 
do not generally appear to be associated with strong electron acceleration, 
perhaps because the Mach numbers of such shocks are much lower than those 
of SNRs and extragalactic radio jets, although Anderson et al. (1979) have 
published data showing that keV electrons are produced in the vicinity of 
the perpendicular bow shock of the Earth. 

While diffusive shock acceleration (Axford et al. 1977; Krymsky 1977; 
Bell 1978; Blandford & Ostriker 1978) provides an efficient means of gener- 
ating highly energetic electrons from an already mildly relativistic threshold, 
and can operate at oblique shocks as well as parallel ones (Kirk & Heav- 
ens 1989), the "injection" or "pre-acceleration" question remains very open: 
by what mechanisms can electrons be accelerated from background (sub- 
relativistic) energies to mildly relativistic energies (Levinson 1996)? In this 
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paper, we investigate one of the possible answers, which has attractive "boot- 
strap" characteristics. Specifically, we suggest that waves are excited by col- 
lective instability of the non-Maxwellian population of ions reflected from 
a perpendicular shock front, and that these waves damp on thermal elec- 
trons, thereby accelerating them. Such a process was first proposed as a 
candidate acceleration mechanism for cosmic ray electrons by Galeev (1984). 
This model was developed further by Galeev et al. (1995), with the added 
ingredient of macroscopic electric fields implied by the need to maintain 
quasi-neutrality in a plasma with an escaping population of electrons. Mc- 
Clements et al. (1997) carried out a primarily analytical study of electron 
acceleration by ion-excited waves at quasi-perpendicular shocks, which was 
necessarily restricted to quantifying linear regimes of wave excitation and 
particle acceleration, in relation to inferred shock parameters. 

Instabilities driven by shock-reflected ions at SNR shocks have also been 
invoked by Papadopoulos (1988) and Cargill & Papadopoulos (1988) as mech- 
anisms for electron heating, rather than electron acceleration. On the basis 
of a simple analytical calculation, Papadopoulos predicted that strong elec- 
tron heating would occur at quasi-perpendicular shocks with "superhigh" 
Mach numbers (specifically, shocks with fast magnetoacoustic Mach numbers 
M F > 30-40) through the combined effects of Buneman (two-stream) and 
ion acoustic instabilities. In this model the Buneman instability, driven by 
the relative streaming of shock-reflected ions and upstream electrons, heats 
the electrons to a temperature T e much greater than the ion temperature Tf. 
in these circumstances, the ion acoustic instability can be driven unstable if 
there is a supersonic streaming between the electrons and either reflected or 
non-reflected ( "background" ) ions. Using a hybrid code, in which ions were 
treated as particles and electrons as a massless fluid, Cargill & Papadopou- 
los (1988) found that the electron heating predicted by Papadopoulos (1988) 
could occur in a self-consistently computed shock structure. However, as 
Cargill & Papadopoulos point out in the last paragraph of their 1988 paper, 
the use of a fluid model for the electrons means that hybrid codes cannot 
be used to investigate electron acceleration. Recently, Bessho & Ohsawa 
(1999) have used a particle-in-cell (PIC) code to investigate acceleration of 
electrons from tens of keV to highly relativistic energies at oblique shocks in 
which the electron gyrofrequency fl e exceeds the electron plasma frequency 
u pe . 

An improved theoretical understanding of electron acceleration at shocks 
is desirable not only for intrinsic interest, but also to enable observations 
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of synchrotron and inverse Compton emission to be related quantitatively 
to shock parameters. However, almost all work on particle acceleration has 
concentrated on ions. There are several reasons for this. First, upstream mo- 
mentum and energy fluxes are dominated by ions, and the shock structure 
problem therefore reduces essentially to that of isotropizing the ion distribu- 
tion. Second, much of our understanding is based on the use of hybrid codes, 
in which electrons are represented as a fluid: such codes cannot provide 
information on electron acceleration. However, the very fact that electron 
dynamics does not appear to be important for shock structure allows us to 
separate the two problems: prescribing the ion parameters using the results 
of hybrid code simulations, we can examine in detail physical processes occur- 
ring on electron timescales. This is the approach followed in this paper. We 
describe the results of a fully nonlinear investigation, carried out by large 
scale numerical simulation using a PIC code and backed up by analytical 
and numerical studies, of the underlying plasma physics mechanisms. We 
consider the case ui pe > Q e , which is qualitatively distinct from the strongly- 
magnetized regime investigated by Bessho & Ohsawa (1999). Our primary 
goal is finding a mechanism capable of producing mildly relativistic electrons: 
once they have attained rigidities comparable to those of shock-heated pro- 
tons, they can undergo resonant scattering, and subsequent acceleration to 
relativistic energies can then proceed via the diffusive shock mechanism. Our 
approach enables us to test earlier predictions of both electron acceleration 
(Galeev 1984; Galeev et al. 1995; McClements et al. 1997) and electron 
heating (Papadopoulos 1988; Papadoulos & Cargill 1988) at very high Mach 
number astrophysical shocks. Simulation results are presented for a range of 
reflected ion speeds in Sect. 2; plasma instabilities occurring in the simula- 
tions, and other processes likely to play a role in electron acceleration and 
heating at SNR shocks, are identified in Sect. 3; and the results of these 
investigations are discussed in Sect. 4. 

2 Particle— in— cell code simulations 

To investigate wave excitation and particle acceleration in the vicinity of a 
perpendicular SNR shock we use an electromagnetic relativistic PIC code 
described by Devine (1995). The particle-in-cell principle (Denavit & Kruer 
1980) relies on self-consistent evolution of electromagnetic fields and macropar- 
ticles in sequential stages. Relativistic electromagnetic PIC codes have been 
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used previously to simulate acceleration processes in astrophysical plasmas 
(e.g. McClements et al. 1993; Bessho & Ohsawa 1999). A distinctive feature 
of the code used in the present study is the fact that the energy density of 
electromagnetic or electrostatic fluctuations can be readily determined as a 
function of frequency u>, wavevector k or time t: this greatly facilitates the 
identification of any wave modes excited in a particular simulation. 

The code has one space dimension (x) and three velocity dimensions 
(v x ,v y ,v z ). To model a plasma containing shock-reflected proton beams, 
we construct a simulation box with 350 grid cells in the x-direction and 
with the local magnetic field B oriented in the y-direction. McClements 
et al. (1997) pointed out that, at any given point in the shock foot, there 
are in fact two proton beams, one propagating away from the shock, the 
other towards it. For simplicity, we assume in our PIC model that the two 
beams propagate with equal speeds Ubi_ in opposite directions perpendicular 
to the magnetic field, and that both background ions and electrons have zero 
net drift: thus, the simulated plasma has zero current. Strictly speaking, 
this is unrealistic, since, in self-consistent models of perpendicular shocks, 
the magnetic field magnitude has a finite gradient along the shock normal 
direction, and a finite current is then required by Ampere's law (Woods 1969). 
We will discuss this approximation in Sect. 4. The frame of reference in each 
simulation is the upstream plasma frame: time evolution in the simulation 
can thus be interpreted as spatial variation in the shock foot, with t = in the 
simulation corresponding to the interface between the undisturbed upstream 
plasma and the foot. The size of the foot Lf oot lies approximately in the 
range (0.3 — 0.7)v s /Qi, where v s is the shock speed and fli is the upstream 
ion gyrofrequency (McClements et al. 1997). Thus, if the simulation is to 
be confined to the foot, the duration of the simulation should be no greater 
than 



where Q e is the electron gyrofrequency (the true proton/electron mass ratio, 
1836, was used in the simulations). The simulations reported here lasted for 
either 70 or 135 electron cyclotron periods 2n/Q e . 

The proton beams were assumed to be initially Maxwellian with ther- 
mal speed 5u± = 3 x 10 5 ms _1 (5u± being defined such that the equivalent 
temperature in energy units is mpSu 2 ^, where m p is the proton mass), and 
a range of perpendicular drift speeds Ub± = 3.25t> e0 , 3.5t> e0 , 5v e0 and 6v e0 , 
where v e0 is the electron thermal speed, defined in the same way as 5u± and 
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initially set equal to 3.75 x 10 6 ms _1 (this corresponds to an electron tem- 
perature T e ~ 9.3 x 10 5 K ~ 80 eV). The value chosen for the total beam 
number density, 0.33n e , is consistent with the highest values of this param- 
eter found in hybrid simulations of quasi-perpendicular shocks with Alvenic 
Mach numbers Ma ranging up to about 60 (Quest 1986). Cargill & Pa- 
padopoulos (1988) used Ma = 50 in their hybrid simulation of an SNR shock 
(it was computationally difficult to simulate shocks with higher Ma). The 
density of each beam rib is, accordingly, one sixth of the electron density n e , 
so that the background proton density rii required by charge balance is 0.67n e 
(the background proton thermal speed Vi was set equal to 1.9 x 10 5 ms _1 ). 
The electron plasma frequency u pe /2n and gyrofrequency Q e /2n in our sim- 
ulations were set equal to 10 5 Hz and 10 4 Hz, respectively, corresponding to 
n e ~ 1.2 x 10 8 m~ 3 and magnetic field B~3.6x 10~ 7 T. The ratio u pe /Vt e is 
typically of order 10 2 or higher in HII regions of the interstellar medium. We 
have chosen a relatively low value of this ratio in order to study and compare 
the effects of instabilities occurring on both the uj~} and fi" 1 timescales. 

The electrons, background protons and each proton beam were repre- 
sented, respectively, by 3200, 800 and 7200 particles per cell. The use of a 
relatively small number of background protons per cell is justified by the fact 
that instabilities driven by the proton beams have much higher frequencies 
than noise fluctuations associated with the background protons: large num- 
bers of electrons and beam protons in each cell ensure a level of noise energy 
well below the wave energy produced by the instabilities. In what follows we 
measure time in electron cyclotron periods, using the notation t = Q e t/27r. 
We also define k = kv e o/Q e (only waves propagating in the x-direction are 
represented), a normalized frequency uj = u/Q e , and r = kv±/Q e , v± being 
electron velocity perpendicular to the magnetic field. 

In every simulation, transfer of energy from beam protons to electrons was 
observed, but the power flux between the two species increased dramatically 
when Ub± was raised from 3.5t> e o to 5v e0 . Figure 1 is a time evolution plot of 
perpendicular kinetic energy S± e = J2j m eV±j/2, where m e is electron mass 
and the summation is over all electrons in the simulation box. Since the total 
electron number is constant, £± e can be regarded as a measure of the effective 
perpendicular electron temperature (although it should be stressed at the 
outset that the electrons do not always have a Maxwellian distribution). 
The energy is normalized to its initial value, which was identical in the four 
simulations. When Ub± = 3.25t> e0 and 3.5v e0 (upper plot) the energy increases 
by approximately an order of magnitude in around 60-100 electron cyclotron 
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periods; when Ub± = 5v e0 and 6t> e0 (lower plot) the energy increases by a 
factor of about 40 within t ~ 15 — 30. The perpendicular energies of the other 
two particle populations, again normalized to the initial electron energy, are 
plotted versus time for the case of Ub± = 6v e o in Fig. 2. In the case of the 
beam protons (upper plot), both bulk motion energy and thermal energy are 
included. During the simulation the beam proton energy drops by less than 
1%, while the background proton energy (lower plot) rises by no more than 
about 10% (in the other simulations the perpendicular energies of the two 
ion species changed by even smaller amounts). In absolute terms the energy 
gained by background protons is very small compared to that lost by beam 
protons, with almost all the energy being transferred to electrons: we will 
demonstrate that the beam protons excite an instability which couples them 
efficiently to electrons. 

In all the cases studied, electrons were energized in the direction per- 
pendicular to the magnetic field. The upper plot in Fig. 3 shows, in more 
detail than Fig. 1, the time evolution of £± e (once again normalized to its 
initial value) in the first 25 electron cyclotron periods of the simulation with 
Uf,± = 6f e o- The lower plot shows the time evolution of (eqE x (x, t) 2 /2), where 
Eo is the permittivity of free space, E x (x,t) is the x-component of the elec- 
tric field, and the brackets () denote a spatial average over the simulation 
box. In general, E x is the dominant field component: since propagation in 
the ^-direction only is represented, it follows that the waves excited are pre- 
dominately electrostatic. Henceforth, the term "electric field" refers to the 
re-component. The field has a single value in each simulation box cell: the 
electrostatic field energy density {e E x (x, t) 2 /2) is calculated by summing 
e E x (x, t) 2 /2 over the box and dividing by the number of cells. The energy 
density plotted in the lower frame of Fig. 3 is normalized to the perpen- 
dicular electron energy density at t — 0. The electron energy grows rapidly 
in two main phases, at t ~ 5 and t ~ 14, and then continues to grow at a 
slower rate. The field energy is greatly enhanced at times when the particle 
kinetic energy is growing rapidly: this suggests strongly that the fields are 
involved in particle acceleration. In the case of the wave energy burst at 
t ~ 5, the field energy grows to a level comparable to the electron kinetic 
energy at that time. The energy of the burst occurring at t ~ 14, on the 
other hand, is much lower than that of the electrons. Figure 4 shows the 
time evolution of £± e and field energy in the simulation with uj,± = 3.25t> e o 
The upper plot shows S± e growing on a timescale comparable to the transit 
time of the simulation box through the shock foot. The lower plot shows 
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that electrostatic field activity is again correlated with electron acceleration. 
Figure 4 resembles the second of the two periods of wave growth in Fig. 3 
(at t ~ 14), in that the wave energy is small compared to the electron kinetic 
energy. 

We now consider the distribution of wave amplitudes in wavenumber 
space. Figure 5 shows the time evolution of this distribution in the sim- 
ulation with Ub± = Gv e o. The grey scale shows the base 10 logarithm of the 
wave amplitude obtained by Fourier transforming in space the electric field of 
one of two counter-propagating waves excited by the ion beams. The start of 
the burst in wave energy in the lower plot of Fig. 3 at t ~ 3 can be identified 
with the burst at k ~ 1.8 in Fig. 5. This reaches an amplitude of 35 Vm -1 , 
generating a harmonic at k ~ 3.6. When the peak amplitude is reached there 
is an increase in wave energy at k < 1. The frequency of this low k noise 
is close to the upper hybrid frequency oj u h = (oJ 2 pe + VL 2 ) 1 / 2 . Its appearance 
correlates with the maximum of the first wave burst at t ~ 5 in the lower 
plot of Fig. 3, and with the strong increase of electron kinetic energy in the 
upper plot, suggesting that it arises from a redistribution of wave energy and 
changes in the electron distribution. After t ~ 8, when the initial wave burst 
has disappeared, a more broadband perturbation is generated at k ~ 1.3, the 
mean k decreasing with time. At t = 14 the wave amplitude peaks at about 
16 Vm" 1 : this is considerably lower than the peak electric field of 35 Vm" 1 
in the first burst, but nevertheless strong enough to generate two harmonics 
(at k ~ 2.6 and k ~ 3.9). 

The corresponding plot for the simulation with Ub± = 3.25t> e o is shown 
in Fig. 6. In this case instability occurs at discrete, regularly-spaced values 
of k. Waves with relatively high k (~ 4) are the first to be driven unstable: 
during the course of the simulation, the instability shifts to lower discrete 
wavenumbers. Broadband noise develops at k < 1, as in Fig. 5, but at a 
later time in the simulation (t ~ 35). This appears to be associated with a 
more gradual evolution of the electron distribution than that which occurs 
in the simulation with Ub± = 6f e o- The difference in temporal behaviour 
between Figs. 5 and 6 will be discussed later in this paper. Figures 5 and 6 
show that in both simulations the plasma eventually stabilizes, on a timescale 
which depends on the beam velocity. 

The dependence of wave amplitude on k and t when Ub± = 5v e o is qual- 
itatively similar to Fig. 5: after an intense burst early in the simulation, 
a wave with slowly-varying amplitude is observed to cascade down in k as 
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time progresses. The growth rate of the first wave burst is 20% higher in 
the simulation with Ub± = 6v £ q than it is in the simulation with u b ± = 5v £ q. 
In the former mentioned above, the peak amplitude of the second 

burst is 16 Vm" 1 , at k — 1.25 and t = 14; the corresponding figures for the 
simulation with Ub± = 5v e o are 12 Vm -1 , k = 1.88 and t = 12. The wave 
amplitude distribution in the simulation with u b ± = 3.5v e0 is similar to Fig. 
6: bursts of wave activity occur at discrete k, with the high k modes being 
driven unstable first. 

In principle, it is also possible to determine the time evolution of wave 
amplitude as a function of Cj and k. However, in order to obtain good fre- 
quency resolution it is necessary to average the amplitude over times longer 
than the electron acceleration timescale. Electrostatic waves in the electron 
cyclotron range propagating perpendicular to the magnetic field include, for 
example, electron Bernstein waves, whose dispersion relation depends on the 
electron distribution. Since this is rapidly evolving, it can be difficult to inter- 
pret observed distributions of wave amplitude in Cj and k. However, we have 
found that the most strongly-growing waves in the simulations invariably 
have Cj ~ kv,b±/v e o: one can thus obtain the frequencies of the high intensity 
modes in Figs. 5 and 6 by multiplying k by Ub±/v e o. By this means, it is 
straightforward to verify that the modes excited early in both simulations 
have Cj ~ 10, and hence cj ~ u pe . 

3 Analysis of simulation results 

Short-lived bursts of narrowband wave activity, correlated with rapid in- 
creases in electron kinetic energy, occur for all four values of Ub±/v e0 con- 
sidered above. These bursts appear throughout the simulations with u b ± = 
3.25f e o and Ub± = 3.5t> e o; in the case of u b ± = 5v e o and Ub± = 6t> e o, they 
appear only at early times. In every case, the instability cascades to longer 
wavelengths in the course of the simulation. In order to compare the simu- 
lation results with those given by linear instability analysis (described in the 
next subsection), we determine growth rates for the first wave that interacts 
significantly with the electrons: in such cases one may assume that the elec- 
trons are still represented by a single Maxwellian velocity distribution with 
thermal speed v e o- The simulations provide the wavenumber k of the un- 
stable wave modes and the electric field amplitude E. The real frequencies 
Cj of the unstable wave modes are assumed to be equal to kub±/v e0 . The 
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normalized growth rate 7/fi e is estimated by fitting an exponential to the 
plot of wave amplitude versus time during the period of most rapid growth 
in each simulation. 

The results of this analysis are shown in Table 1. The symbol E m denotes 
the maximum value of E during each simulation. In three of the four simula- 
tions there is a period of wave growth which can be described accurately as 
exponential. In each case, the growth rate falls to zero, and the wave decays: 
an example of this behaviour, for the case of Ub± = 6t> e o, is shown in Fig. 
7, where wave amplitude (defined as in Figs. 5 and 6) at k — 1.8 is plot- 
ted versus normalized time. A possible reason for wave collapse (observed 
in all four simulations) will be discussed later in this paper. In the case of 
Ub± = 3.25t> e0 , the mode referred to in Table 1 {k ~ 3.6) is the second to 
be destabilized in that simulation. It appears to grow linearly rather than 
exponentially: for this reason, no figure is given for its growth rate. The 
first mode to be destabilized in this simulation, with k ~ 3.9, does not grow 
to a large amplitude (compared to the noise level), and so it is difficult to 
determine its growth rate. Later in the paper we will present evidence of 
wave-wave interaction between the second mode excited (k ~ 3.6) and the 
third mode excited (k ~ 3.3), which may help to explain the linear growth 
of the latter. 

Table 1. Parameters of highest intensity wave mode in each simulation. 



Ub_l_/v e0 


k 




7M 


E m (Vm- 1 ) 


6.0 


1.8 


10.8 


0.24 


35 


5.0 


2.15 


10.7 


0.2 


23 


3.5 


3.3 


11.6 


0.05 


2.5 


3.25 


3.6 


11.7 




1.6 



Let us now examine whether the growth rates derived from the PIC simu- 
lations in Table 1 and Fig. 7 can be reproduced using linear stability theory. 

3.1 Linear stability analysis 

The appropriate dispersion relation for electrostatic, perpendicular-propagating 
waves with frequencies in the electron cyclotron range and above, excited by 
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an ion beam with a Maxwellian distribution in v±, is (Melrose 1986) 




ft [1 + ( b Z(( b )} 
k 2 Su 2 , 



pe 



,2 -\ e oo 



uj — iVl e 



Pie 



0, 



(2) 



oo 



where: uj p i, uj p b are the background and beam ion plasma frequencies; Z is 
the plasma dispersion function, with argument = (oj — kub±)/kSu±; and 
Ie is the modified Bessel function of the first kind of order I with argument 
A e = T e k 2 /(m e Ql). Both species of ion, having a much longer cyclotron 
period than the electrons, can be treated as unmagnetized particles on the 
timescales of interest here. Strictly speaking, there should be a term in Eq. 
(2) for each of the two proton beams, but since they have mean perpendicular 
speeds of opposite sign, and uj ~ kub± is a prerequisite for wave-particle 
interaction, we need only consider one of them. 

Solutions of Eq. (2) for complex uj in terms of real k can be readily 
obtained numerically, and compared with the simulation results in Table 1. 
In Figs. 8 and 97 = Im(d>) is plotted versus k for parameters corresponding 
to the initial conditions of the simulations with Ub± = 6v e0 and Ub± = 3.25t> e0 . 
In the former case it can be seen that strong instability drive occurs at 
k ~ 1.8 with maximum growth rate 7 ~ 0.25, as observed early in the 
simulation (Fig. 5 and Table 1). The unstable real frequencies range from 
uj ~ 8 to uj ~ 10.8, and are thus clustered around the dimensionless electron 
plasma frequency {uj = 10). The main instability appears to be essentially 
unaffected by cyclotronic effects: the growth rate does not depend on how 
close the frequency is to a cyclotron harmonic. There are, however, two much 
weaker instabilities at k < 1, which are narrowband in both k and uj, the real 
frequencies lying just below the second and third cyclotron harmonics. In 
Fig. 9 instability occurs at k ~ 3—4, the corresponding real frequencies again 
clustering around uj pe . In this case, however, the instability is modulated by 
cyclotronic effects, as in the simulation. Instability again occurs at k < 1, 
with real frequency uj ~ 1.8. 

The mode appearing early in the simulation with Ub± = Qv c q arises from 
a two-stream instability (Buneman 1958). This can be driven by ions drift- 
ing relative to electrons in an unmagnetized plasma: it can also occur in a 
magnetized plasma, with ions drifting across the field, if uj pe /Q e is sufficiently 
large, and the instability drive is sufficiently strong. Electrons as well as ions 
are then effectively unmagnetized and the appropriate dispersion relation is 
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(Melrose 1986) 



CO 



2 2uj* b {l + CZ(( b )] 2^ 2 e [l + C^(Ce 



-~pi pb L 1 ^~V^"/J — pe L ^ _ r, /o\ 



where ( e = u/kv e Q. In the frequency regime of interest here (a; ~ u; pe ), it can 
be shown easily that the background ion term in Eq. (3) can be neglected. 
Letting the thermal speeds of the two remaining species tend to zero, Eq. 
(3) reduces to 



1 LP£ 1P£ = o . (4) 

This differs slightly from the original two-stream dispersion relation analysed 
by Buneman (1958) in that the ions, rather than the electrons, have a finite 
drift speed. Buneman's equation becomes identical to Eq. (4) under the 
transformation u — > ku b j_ — oj: using this, we can infer from Buneman's 
analysis that Eq. (4) has a root uo = ujq + ry, where real frequency uq and 
growth rate 7 are given approximately by 

luo ~ ku b± - cuTiuli 3 cos 4/3 9, (5) 



J pb w pe 

^ ~oj 2 p / b 3 uj 1 p / e 3 cos 1/3 9sm9, (6) 

9 being a parameter whose value depends on (ku b ^— uu pe ) / Jt^ uu^ (Buneman 
1958): it is straightforward to verify that the strongest wave growth occurs 
when 9 = n/3, which corresponds to ku b ± = oo pe . If ku b ± ^> J^uj p {? cos 4 / 3 9, 
it follows from Eq. (5) that uo ~ ku bl _ and so the strongest drive occurs at 
uj ~ uj pe . However, since 9 can have a range of values, the instability has 
finite bandwidth, extending to frequencies significantly below u pe . Solving 
the full Buneman dispersion relation [Eq. (4)] with u b ± = Qv e o, we obtain 
results which are almost identical to those obtained from the magnetized 
dispersion relation [Eq. (2)], except, of course, that the cyclotronic features 
at k < 1 in Fig. 8 do not appear. Even in the case of u b ± = 3.25f e0 , 
Eq. (4) yields instability at about the same wavenumbers and frequencies as 
Eq. (2) [although the growth rates are somewhat higher in the case of Eq. 
(4)]. The essential difference between Figs. 8 and 9 is that the lower beam 
speed in the latter yields lower growth rates: when 7 is sufficiently small, the 
gyromotion of an electron in one wave growth period cannot be neglected, 
and the instability is modified by cyclotronic effects. However, the instability 
remains Buneman-like in character. 
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Further analysis of Eq. (2) indicates that the instability growth rate 
is a slowly-decreasing function of 5uj_/ubj_: in the case of u b ± = Qv 6 q, for 
example, the maximum growth rate is around 0.06f2 e when Su±/ub± = 0.3. 
The Buneman instability is thus robust, in the sense that its occurrence 
is not critically dependent on the velocity-space width of the reflected ion 
distribution. In any event, the values of 5u^_/ubi_ used in our PIC simulations 
are broadly consistent with reflected beam ion distributions occurring in the 
hybrid simulations of Cargill & Papadopoulos (1988). 

The instabilities at k < 1 in Figs. 8 and 9 arise from the interaction of 
a beam mode (uj ~ kubi_) with electron Bernstein modes. The existence of 
such instabilities can be inferred analytically by taking the limit of Eq. (2) 
for cold beam and background protons: 

i .2 , ,2 , ,2 „-X e oo o2 T 

uj 2 (uj - ku b± ) 2 uj \ e uj - £VL e 

In the absence of the proton beam term, Bernstein mode solutions of Eq. (7) 
have frequencies which approach £fl e (£ = 1,2, 3, ...) as k — > oo and (£+l)Q e 
as k — > (the long wavelength limit is different for frequencies equal to 
or greater than the upper hybrid frequency u u h, which in the case of the 
simulations presented in Sect. 2 is about 10Q e ) . When rib 7^ 0, approximate 
analytical solutions of Eq. (7) can be obtained by setting u = kub± + «7 and 
solving perturbatively for 7 in certain limits. For example, letting A e — > 
and assuming that uj does not lie close to a harmonic of Q e , we obtain 



J_ 




(8) 



Instability, corresponding to real 7, thus requires uj > Q e . For uj sufficiently 
close to £Q e , the electron term in Eq. (7) is dominated by the £-th harmonic, 
and instead of Eq. (8) we obtain 



tt e tt e \miJ \n e J (£I e ) 1/2 



UJ 



(9) 



Numerical solutions of Eq. (2) for uj ~ f2 e are broadly consistent with Eqs. 
(8) and (9). In both these cases the growth rate scales as {m e lm, p ) l l 2 : in 
contrast, the Buneman growth rate [Eq. (6)] scales as (m e /m p ) 1 ^ . This 
helps to explain the fact that in Figs. 8 and 9 the Buneman instability 
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is stronger than the lower frequency Bernstein instability. It should also 
be noted that most astrophysical plasmas have a higher ratio of electron 
plasma frequency to gyrofrequency that that assumed in the simulations 
(up e /Q e = 10). Normalized to Q e , the Buneman growth rate scales as u pe /fl e , 
and so the instability is less likely to be modified by cyclotronic effects when 
uo pe /VL e > 10. The electron Bernstein modes exist because T e is finite: thus, 
the transition from Buneman to electron Bernstein instability depends on the 
value of v e Q. If the initial electron temperature in the simulations had been 
lower than 80 eV, the Buneman instability would, again, have been affected 
to a lesser extent by cyclotronic effects. 

3.2 Nonlinear effects 

Figure 1 shows strong increases in electron kinetic energy perpendicular to 
the magnetic field in all four simulations. There is a strong correlation be- 
tween acceleration and wave excitation via the Buneman instability (Figs. 
3 and 4). Although such waves can energize electrons via Landau damping 
(Papadopoulos 1988), one would expect this process to be of limited effective- 
ness when, as in the present case, the waves are propagating perpendicular to 
the magnetic field and have a growth rate which is comparable to or less than 
f2 e . It is likely therefore that the very strong electron acceleration observed 
in the simulations is due at least in part to nonlinear processes. 

As noted previously, the second mode to be excited in the simulation 
with Ub± = 3.25t> e o does not undergo an exponential growth phase. Fig- 
ure 10 shows the time evolving wave amplitudes of this mode, at k — 3.6 
(upper plot), and the third mode to be excited, at k — 3.3 (lower plot). 
The amplitude at k — 3.6 grows linearly up to t ~ 27, and then collapses. 
The amplitude at k — 3.3 grows exponentially from t ~ 7 to t ~ 14, with 
■y/Q e ~ 0.04: this is close to the growth rate at k — 3.3 found by linear 
stability analysis (Fig. 9). The amplitude remains constant until t ~ 27, and 
then grows linearly until t ~ 35. The linear growth of the wave at k — 3.3 
thus correlates strongly with the collapse of the wave at k — 3.6: this sug- 
gests wave-wave coupling. The linear growth of the wave at k — 3.6 may, in 
turn, be correlated with the decay of the first mode to be excited, at k ~ 3.9 
(see Fig. 6): the latter has the highest growth rate of waves in this range, 
according to linear stability analysis (Fig. 9). 

We now consider specific nonlinear processes which may be occurring in 
the simulations. Karney (1978) examined the nonlinear interaction of large 
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amplitude electrostatic lower hybrid waves with ions. The particle motion is 
described by a normalized Hamiltonian 



1 1 

h= -{p x + y) 2 + -p y 2 - a sin (y - [it), (10) 

where: y, z are, respectively, the wave propagation and magnetic field direc- 
tions; the canonical momentum components p x , p y are normalized to mfl/k, 
where m and Q are particle mass and gyrofrequency; the position variable y 
is normalized to 1/k; /i is wave frequency in units of Q; and a is given by 

E being the wave electric field amplitude. The first two terms in the Hamilto- 
nian describe the motion of the particle in the magnetic field; the third term 
arises from the electrostatic wave. The system can be regarded as consist- 
ing of two harmonic oscillators: one associated with the particle gyromotion 
around B, the other with the wave. These oscillators are coupled by the 
parameter a, the value of which determines the extent to which the system 
phase space is regular or stochastic. Karney (1978) solved the Hamiltonian 
equations corresponding to Eq. (10) for a range of initial conditions, plot- 
ting normalized Larmor radius r = kv±/Q versus wave phase angle at the 
particle's position, for successive transits of the particle through a particu- 
lar gyrophase angle. Particle trajectories were thus represented as discrete 
sets of points in (r, 0) space. For sufficiently small values of a, all particles 
have regular orbits, represented by smooth curves in (r, 0) space, spanning 
all values of and with only small variations in r. When a exceeds a certain 
threshold a , "islands" appear in (r, 0) space, within which particle trajec- 
tories are confined. Further increases in a lead to the formation of more 
islands, which cause the phase space to become stochastic: at sufficiently 
large a > a c , the system phase space is completely stochastic, with no reg- 
ular orbits remaining. The initial electron distributions in our simulations 
decrease monotonically in v±: in such cases stochasticity in phase space tends 
to favour particle diffusion to larger velocities, i.e. acceleration. 
Karney (1978) obtained the following analytical estimate for a : 



a 



r(u/n - i) 



£(d/dr)J e (r) 



(12) 
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where £Q is the cyclotron harmonic closest to uo and Ji is the Bessel function 
of order I. Karney's analysis does not explicitly involve a particular type of 
wave or particle, or a specific particle distribution function. The results can 
thus be applied to the case of electrons interacting with electrostatic waves 
excited by the Buneman instability, in which case m = m e and Q = Q e . 
Combining Eqs. (11) and (12) and using the identity 

^-Mr) = -J e (r)-J £+1 (r), (13) 
ar r 

we infer that the critical electric field E = Ei for island formation in (r, <fi) 
space is 

v_ 1 B^-i\ 

1 t\lJ e (r) - J t+1 (r)\- U4j 

In general, it is not possible to determine analytically an expression for the 
electric field amplitude E = E c corresponding to a = a c , above which the 
phase space becomes completely stochastic. Karney obtained an empirical 
expression for a c , based on numerical calculations with particular values of \x 
and r, which may not be applicable to the simulation results discussed here. 
However, island formation is a first step in the destruction of regularity in 
the system phase space, and Ei can be regarded as an approximate threshold 
for stochasticity: electric field amplitudes which are significantly higher than 
E; L will convert regular orbits at a particular v± into stochastic ones. 

In the cases Ub± = 5v e0 and 6v e0 , linear stability analysis indicates that 
wave growth occurs across a range of frequencies uj ~ u pe , which includes 
cyclotron harmonics: in such cases fi = £, and any non-zero wave amplitude 
E will cause islands to be formed. In the case of the lower two beam speeds, 
the unstable frequencies lie between cyclotron harmonics, and Ei is thus 
always finite. Table 2 lists the values of Ei derived from Eq. (14) that are 
required for comparison with the highest intensity wave mode excited in each 
simulation. The actual peak electric fields of these waves are given in Table 
1. 
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Table 2. Values calculated for E { using the wave parameters given in Table 
1. 



v±/v e0 


closest £ 




kv±/Q e 


Ei (Vm- 1 ) 


6.0 


11 


0.2 


10.8 


1.8 


5.0 


11 


0.3 


10.7 


2.3 


3.5 


12 


0.4 


11.6 


2.1 


3.25 


12 


0.3 


11.7 


1.4 



Comparing Tables 1 and 2, we see that waves are excited with amplitudes 
exceeding E i in all four cases. For the simulation with u^± = 3.25t> e0 , E/Ei ~ 
1.1. This ratio rises to 1.2 for Ub± = 3.5t> e0 , 10 for Ub± = 5v e0 , and 19 for 
Ub± = 6f e0 . In the latter two cases, as we have seen, waves are excited with 
uj = £Q e , for which island formation occurs regardless of the value of E. The 
fact that E m /Ei ^> 1 at higher values of Uf,± indicates that the phase space in 
these simulations is characterized by strong stochasticity. The waves rapidly 
collapse, however, soon after the onset of strong electron acceleration. In 
the other two simulations, the peak amplitudes are only just sufficient for 
island formation to occur, and it is likely that little stochasticity occurs in 
the system phase space. The waves excited in these simulations decay more 
gradually than those produced at higher Uf,±. 

We now consider possible explanations for two of the results noted above: 
the sharp rise in wave amplitude when the beam speed is raised from 3.5v £ q 
to 5v e o; and wave collapse, which occurs in all four simulations but is particu- 
larly rapid in the two simulations with higher uj,±. As far as the dependence 
of wave amplitude on Ub± is concerned, the first point to note is that the 
unstable waves all satisfy uj ~ Uf,±k. In each simulation the total number 
of computational particles is, of course, finite, the Maxwellian electron ve- 
locity distribution being initialized up to v± ~ 5v c q. Thus, the beams with 
Ub± = 5v £ q and 6t> e o excite waves with phase velocities exceeding the veloc- 
ity of any electron in the simulation: this is not so in the simulations with 
Ub± = 3.25t> e0 and 3.5v e0 . The minimum electron velocity required for wave- 
particle interactions is the phase velocity of the wave: thus, only the slow 
beams can excite waves capable of interacting with electrons at the start of 
the simulations. The wave-particle interaction results in electron accelera- 
tion, the energy for this being drawn from the wave. This energy loss may 
account for the relatively low peak electric field amplitudes of waves excited 
by the slow beams. 
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The waves generated by the fast beams, on the other hand, cannot ini- 
tially interact resonantly with the electron population, and so their ampli- 
tudes can grow to levels much higher than Ei. Sufficiently high wave am- 
plitudes can activate a second acceleration mechanism, which arises from 
particle trapping in the wave electric field (Karney 1978): electrons with 
an initially monotonic decreasing distribution are re-distributed uniformly 
within the trap, the result being a net increase in kinetic energy. The wave 
can trap electrons with perpendicular velocities differing from the wave's 
phase velocity by up to v tr , where 



For uj,± = 5t> e o, the maximum electric field is 23 Vm and the wavenumber 
k is 5.7 x 1CT 3 x 27rm~ 1 . In this case v tT = 1.1 x 10 7 ms~ 1 ~ 2.8t> e0 . For 
u b ± = 6f e0 , the maximum field is 35 Vm -1 and fc = 4.8x 1CT 3 x 27rm _1 , so 
that v tr = 1.4 x 10 7 ms _1 ~ 3.8t> e0 . The waves excited in the simulations with 
higher beam speeds can thus trap electrons deep within the electron ther- 
mal population: a large number of electrons can then be pre-accelerated to 
velocities comparable to the wave's phase velocity, with further acceleration 
taking place via the stochastic mechanism discussed previously. A two-stage 
process of this type was proposed by Karney (1978). Whereas the first burst 
of wave activity in Fig. 3 contained more energy than the electron popula- 
tion, the energy in the second burst was much lower than the perpendicular 
electron kinetic energy by that stage of the simulation. This may have been 
due to the first burst resulting in trapped electrons populating the region of 
phase space at v± ~ Ubj_, via the trapping mechanism. The perpendicular 
electron velocity distribution would then be considerably broader than it was 
initially, with an effective thermal speed v e > v £ q. The beam distribution, on 
the other hand, did not change significantly during the simulation (see Fig. 
2), and so Uf,±/v e < Ubj_/v e o. The situation would then be similar to that of 
the simulations with lower beam speeds, in which electrons can immediately 
absorb energy from waves with u ~ kub±, and one would expect any subse- 
quent wave burst to have a peak energy much lower than that of the electron 
population, as observed. 

With regard to the second observation, wave collapse, it is interesting to 
note that in every case the wave amplitude falls to a level well below E^. 
intuitively, one would have expected the waves to cease interacting with elec- 
trons, and hence to reach a steady-state level, when their amplitudes had 




(15) 
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fallen below E { . The collapse may be associated with changes in the disper- 
sion characteristics of the wave mode resulting from strong particle accelera- 
tion. Karney (1978) justified his Hamiltonian approach by considering only 
stochastic regions of phase space, at particle speeds (and hence wave phase 
speeds) greatly exceeding v e0 . The stochastic regions thus lie in the high ve- 
locity tail of the initial Maxwellian electron distribution, and most electrons 
are not initially affected by the wave-particle interaction. However, in the 
simulations with Ub± = 5t> e o, Qv e o the reduction in Ub±/v e noted above means 
that perpendicular electron speeds are no longer small compared to the wave 
phase speed, and we find that there is a transition from the pure Buneman 
instability shown in Fig. 8 to the more complicated instability shown in 
Fig. 9: the latter, as we have discussed, has a Buneman-like envelope, but 
also has cyclotronic features, and in fact linear stability analysis shows that 
the variation of uo with k in this case is characteristic of the beam/electron 
Bernstein mode discussed in Subsect. 3.1. As Ub±/v e falls, the maximum 
growth rate drops considerably, but remains positive if the electrons retain 
a Maxwellian distribution. However, as we now demonstrate, the electron 
distributions occurring in the simulations are often far from Maxwellian. 

3.3 Particle distributions 

From the simulation results we have evaluated the distribution of perpendic- 
ular electron speeds f(v±), defined such that 



where N e is the total number of electrons in the simulation. With this defini- 
tion, a Maxwellian velocity distribution is of the form v±e~ v ^^ 2Ve , decreasing 
monotonically for v± > v e . One advantage of plotting a distribution in this 
way is that the thermal speed of a Maxwellian can be readily identified graph- 
ically, being the speed at which df/dv± = 0. In Fig. 11 f(v±) is plotted for 
t = 0, 45, 90 and 135 in the simulations with Ub± = 3.25t> e0 (dash-dotted 
curves) and Ub± = 3.5t> e o (solid curves). The two curves are identical for 
t = 0, since the same initial electron temperature is used in all four sim- 
ulations. At t = 45 the proton beams have generated hot electron tails, 
peaking at v± ~ 4t> e0 (ub± = 3.25v e o) and v± ~ 6v e o (ub± = 3.5t> e o). The 
maximum electron speeds in the two cases are 7v £ q (ub± = 3.25t> e o) and 10v £ q 
{ub±_ = 3.5t> e0 ). At t = 90 the slower beam has produced a local maximum 




(16) 
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i n f( v ±) a ^ 5f e0 , and a high velocity cutoff at 10v e0 . The local maximum 
has become less pronounced at t — 135. In the case of Ubj_ = 3.5t> e o a lo- 
cal maximum can be seen at t — 45: by t = 90, however, the distribution 
is monotonic decreasing above a speed only slighly higher than the initial 
thermal speed. By the end of this simulation f(v±) extends up to 12t> e0 . 

In Fig. 12 f(v±) is shown for t — 0, 20, 40 and 70 in the simulations with 
Ub± = 5t> e o (dash-dotted curves) and Ub± = Qv £ q (solid curves). At t = 20 
the two distributions have local maxima at v± ^> v e o, as in the second frame 
of Fig. 11: for Ub± = 5v e o the distribution peaks locally at v± ~ 10t> e o and 
falls to zero at v± ~ 18t> e0 . The corresponding figures at the same stage of 
the simulation with Ub± = 6v e0 are 12v e0 and 22v e0 . By t = 40, the local 
maxima still exist and, indeed, the bumps-on-tail containing these maxima 
actually comprise most of the electron population in both cases. By this 
time the high velocity tails extend to v ± ~ 25 — 30t> e o- Local maxima close 
to the original thermal speeds v e o still exist, but these have disappeared by 
t = 70. The strong wave bursts in these simulations occur before t — 20: after 
this time, a weaker, more broadband instability occurs at lower k, but still 
with Co ~ kub±/v e0 . To model this instability, we can approximate the solid 
curve at t — 20 in Fig. 12 by superposing two Maxwellians, with thermal 
velocities v ec = v e o, v e h = 10v e o and densities n ec = 0.38n e , n e h = 0.62n e . 
The solid curve in Fig. 13 shows the true distribution at t — 20; the dashed 
curve shows the bi-Maxwellian fit to this distribution. The match is not 
exact, but is close enough to suggest that we can model wave excitation at 
this stage of the simulation by solving a modified version of the dispersion 
relation [Eq. (2)], with the parameters of the dashed curve defining the 
electron distribution. Results indicate an electron Bernstein instability with 
maximum growth rate 7 ~ 0.08fi e at k ~ 1.4, Co ~ 8.2f2 e : the wavenumbers 
are consistent with those of fluctuations appearing at t > 20 in Fig. 5. 

The electron distributions at t = 70 in the simulations with Ub± = 5v e o 
and Ub± = 6f e o (Fig. 12) can both be approximated by single Maxwellians, 
respectively with v e ~ 8t> e o and v e ~ 12f e0 . The proton beam-excited Bune- 
man instability can thus produce electron distributions whose perpendicular 
thermal speeds exceed the velocities of the ion beams which produced them. 
The fastest-moving electrons have v± ^> Ub±. This phenomenon, observed in 
all four simulations, is further strong evidence for nonlinear processes playing 
a role in electron acceleration: in the quasi-linear limit, unmagnetized elec- 
trons of a particular speed v can only interact with waves whose phase speed 
equals v , and the range of wave phase speeds is determined in turn by the 
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ion beam speeds. In the case of Ub± = Qv £ q, the final electron temperature is 
about 11.5 keV: this is easily sufficient to account for thermal X-ray emission 
observed from SNRs such as Cas A (Papadopoulos 1988). Individual electron 
energies up to several tens of keV were obtained in the simulations. 

4 Conclusions and Discussion 

Using particle-in-cell (PIC) simulations and linear stability theory, we have 
shown that electrostatic waves in the electron plasma range, excited by ions 
reflected from a high Mach number perpendicular shock, can effectively chan- 
nel the free energy of the shock into electrons. Such shocks are known to be 
associated with SNRs, and the processes investigated in this paper may thus 
help to account for both X-ray thermal bremsstrahlung and the creation of 
"seed" electron populations for diffusive shock acceleration to MeV and GeV 
energies in such objects. The simulation results provide confirmation of a 
proposal by Papadopoulos (1978) and Cargill & Papadopoulos (1978) that 
streaming between reflected ions and upstream electrons can give rise to a 
strong Buneman instability. Whereas these authors assumed that the sole 
effect of the Buneman instability would be electron heating, the PIC sim- 
ulations show acceleration - the development of strongly non-Maxwellian, 
anisotropic features in electron velocity distributions. The maximum elec- 
tron velocities are considerably higher than those expected on the basis of 
quasi-linear theory: this implies that nonlinear wave-particle interactions 
are contributing to the acceleration. When the beam speed is greater than 
about four times the initial electron thermal speed, thermalization of the 
electron population is observed after saturation of the Buneman instability, 
the final electron temperature being of the order of 10-12 keV. 

It is possible that the acceleration process identified in this paper may 
be relevant to oblique shocks as well as strictly perpendicular ones. A 
necessary requirement is the presence of reflected ion beams, which have 
been observed (Sckopke et al. 1983) upstream of both quasi-parallel and 
quasi-perpendicular regions of the Earth's bow shock (the term "quasi- 
perpendicular" is conventionally used to describe a shock at which the angle 
9bu between the shock normal and the upstream magnetic field is greater 
than 45°). Leroy et al. (1982) used a hybrid code to simulate ion reflec- 
tion at shocks with a range of values of finding very similar results 
for 6bu — 80° and 0Bn = 90°. They inferred from this that hybrid sim- 
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ulations which use strictly perpendicular geometry may be compared with 
observational data of quasi-perpendicular shocks. Additional necessary re- 
quirements for electron acceleration via the Buneman instability are that the 
projection of the reflected ion beam velocity onto the plane perpendicular to 
the upstream field exceed the upstream electron thermal speed, and that the 
plasma be weakly magnetized, in the sense that u pe > Q e . When these con- 
ditions are satisfied locally, the Buneman instability will occur. Whether this 
instability is sufficient to produce significant electron energization at oblique 
shocks as well as perpendicular and nearly perpendicular ones remains to be 
demonstrated, however. In the simulations presented in this paper, accel- 
eration occurred on timescales shorter than an ion cyclotron period. It was 
not necessary to represent the shock foot structure, since this has dimensions 
of the order of a reflected ion Larmor radius, and for this reason 9bu is not 
explicitly a parameter in our model. Leroy et al. (1982) have noted, how- 
ever, that extrapolations of results obtained for nearly perpendicular shocks 
(80° < 9 Bn < 90°) to more oblique shocks should be treated with caution, 
since the physical processes governing the shock structure can be expected 
to change as is reduced. 

The simulation results and our analysis of them suggest that the damp- 
ing of waves in the electron plasma range at perpendicular SNR shocks could 
provide a solution to the cosmic ray electron injection problem. Although 
wave-particle interactions at such shocks have been invoked previously in 
this context (Galeev 1984; Papadopoulos 1988; Cargill & Papadopoulos 1988; 
Galeev et al. 1995; McClements et al. 1997), the simulation results presented 
here contain several new features. These include: the acceleration, rather 
than heating, of electrons via the Buneman instability; the acceleration of 
electrons to speeds exceeding those of the shock-reflected ions producing the 
instability; and strong acceleration of electrons perpendicular to the mag- 
netic field. The wave-particle mechanisms proposed by Galeev (1984) and 
McClements et al. (1997) gave rise to electron acceleration primarily along 
the magnetic field. Diffusive shock acceleration, which is probably essential 
for the production of ultra-relativistic electrons, can occur when the elec- 
trons have magnetic rigidities comparable to those of ions flowing into the 
shock. Since, however, the diffusive shock mechanism requires electrons to 
be rapidly scattered, its efficacy does not depend sensitively on the initial 
pitch angle distribution. The geometry of the simulations described here 
(B perpendicular to the one space dimension) excludes the possibility of ac- 
celeration by electrostatic waves along the parallel direction. The present 
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model is complementary to those of Galeev (1984), Galeev et al. (1995) and 
McClements et al. (1997), in that it provides an alternative means of ener- 
gizing electrons at perpendicular shocks. At a real SNR shock, perpendicular 
acceleration via the Buneman instability and parallel acceleration via wave 
excitation at u < Q e are both likely to occur. PIC simulations in two space 
dimensions would make it possible to assess quantitatively the relative con- 
tributions of these two types of instability to electron energization under a 
range of conditions. 

We discuss finally our neglect of the finite plasma current present in the 
foot region of perpendicular shocks. This approximation does not appear 
to have introduced unrealistic elements into our simulation results, except 
insofar as the absence of a finite drift between the electrons and background 
protons removes a possible source of drive for the ion acoustic instability, in- 
voked as one of the electron heating mechanisms in the model of Papadopou- 
los (1988). However, if the background protons and electrons flowing into the 
shock foot are approximately isothermal, it seems unlikely that any instabil- 
ity excited by their relative streaming could result in a significant transfer of 
energy from one species to the other. Another possibility is that ion acous- 
tic instability could result from the streaming of beam protons relative to 
electrons. This would require, however, the electron temperature to be ex- 
tremely large compared to the beam proton temperature (ion Landau damp- 
ing strongly suppresses the instability when the temperatures are equal): 
even in the simulation which produced an effective electron temperature of 
80 times the initial temperature, this condition was not satisfied. Thus, the 
simulation parameters were such that the ion acoustic instability could not 
occur. It is interesting to note, however, that the curve corresponding to 
Ub± = 6v e o in the lower frame of Fig. 1 bears a certain resemblance to a 
curve in Fig. 5 of Papadopoulos' 1988 paper (Fig. 5), showing schematically 
the predicted variation of electron temperature with distance in a quasi- 
perpendicular shock foot: in both cases, there is a very rapid rise in the 
total electron energy, resulting from strong Buneman instability, followed by 
a more gradual rise, which coincides in the simulations with the excitation of 
electron Bernstein modes. The latter may play a role in the simulations which 
is similar to that of the ion acoustic mode in the model of Papadopoulos. 
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Figure 1: Total electron perpendicular kinetic 
energy, normalized to its initial value, versus 
simulation time in electron cyclotron periods 
27r/f2 e , for several values of Ub±/v e o- Energy 
transfer to electrons is much more rapid at Ub± = 
5v e0 and u b ± = 6v e0 than it is at lower beam 
speeds. 
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Figure 2: Normalized beam proton (upper plot) 
and background proton (lower plot) perpendic- 
ular kinetic energies (both normalized to ini- 
tial electron thermal energy) versus simulation 
time for the simulation with Ub± = 6t> e o- The 
beam proton energy drops by less than 1%; the 
background proton energy increases by approx- 
imately 10%. 
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Figure 3: Time evolution of perpendicular elec- 
tron kinetic energy (upper plot) and electro- 
static field energy (lower plot) in the simulation 
with Ub± = 6v e0 . 



27 



60 80 100 120 

Normalized time 



60 80 100 120 

Normalized time 



Figure 4: Time evolution of perpendicular elec- 
tron kinetic energy (upper plot) and electro- 
static field energy (lower plot) in the simulation 
with Ub± = 3.25f e0 . 
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Figure 5: Base 10 logarithm of electric field am- 
plitude (Vm _1 ) versus k and t in the simulation 
with Ubi_ = 6i>eo- 
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Figure 6: Base 10 logarithm of electric field am- 
plitude (Vm _1 ) versus k and t in the simulation 
with Ubi_ = 3.25f e o. 
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Figure 7: Time evolution of electric field ampli- 
tude (Vm -1 ) at k = 1.8 after the onset of in- 
stability in the simulation with Uf,± = 6v e0 . The 
horizontal line is Ei for this wave; the diagonal 
line shows an exponential fit to the growth rate 
with 7/fi e = 0.24. 
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Figure 8: Predicted linear growth rates of waves 
with u > Q e when the beam speed is 6v c q. The 
other dispersion relation parameters correspond 
to the initial conditions of all four simulations. 
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Figure 9: Predicted linear growth rates of waves 
with uj > Q e when the beam speed is 3.25f e o- 
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Figure 10: Time evolution of electric field am- 
plitude (Vm _1 ) at k — 3.6 (upper plot) and 
k = 3.3 (lower plot) after the onset of insta- 
bility in the simulation with Ub± = 3.25t> e o- The 
vertical lines indicate a period in which there is 
an anti-correlation between the two wave am- 
plitudes. 
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Figure 11: Normalized perpendicular electron 
speed distributions at various times in the sim- 
ulations with = 3.25t> e0 (dash-dotted lines) 
and Ub± = 3.5f e0 (solid lines). 
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Figure 12: Normalized perpendicular electron 
speed distributions at various times in the simu- 
lations with Ubi_ = 5t> e o (dash-dotted lines) and 
Ub± = 6f e o (solid lines). 
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Figure 13: Normalized perpendicular electron 
speed distribution at t = 20 in the simulation 
with = 6v e0 . The solid curve is the distribu- 
tion sampled in the simulation; the dashed curve 
shows a bi-Maxwellian fit. 
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